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Abstract 

We consider a branching particle system where particles reproduce according to the pure 
birth Yule process with the birth rate A , conditioned on the observed number of particles to 
be equal n. Particles are assumed to move independently on the real line according to the 
Brownian motion with the local variance a 2 . In this paper we treat n particles as a sample 
of related species. The spatial Brownian motion of a particle describes the development of 
a trait value of interest (e.g. log-body-size). We propose an unbiased estimator R 2 of the 
evolutionary rate p 2 = a 2 /X. The estimator R 2 is proportional to the sample variance S 2 
computed from n trait values. We find an approximate formula for the standard error of R 2 
based on a neat asymptotic relation for the variance of S 2 . 

(Keywords: Branching Brownian motion, conditioned branching process, tree-free phyloge- 
netic comparative method, quantitative trait evolution, Yule process) 

1 Introduction 

Biodiversity within a group of n related species could be quantified by comparing suitable trait 
values. For some key trait values like log body size, researchers apply the Brownian motion 

model proposed by Felsenstein [1985]. It is assumed that the current trait values (Xj , . . . ,X„ ) 
have evolved from the common ancestral state Xq as a branching Brownian motion with the local 
variance a 2 . Given a phylogenetic tree describing the ancestral history of the group of species 
the Brownian trajectories of the trait values for sister species are assumed to evolve indepen- 
dently after the ancestor species splits in two daughter species. The resulting phylogenetic sam- 
ple (X , . . . ,xjfi) consists of identically distributed normal random variables with a dependence 
structure caused by the underlying phylogenetic signal. 

A mathematically appealing and biologically motivated version of the phylogenetic sample 

model assumes that the phylogenetic tree behind the normally distributed trait values (X^ , . . . , X„ ) 
is unknown. As a natural first choice to model the unknown species tree, we use the Yule process 
with birth rate A [see Yule, 1924]. Since the phylogenetic sample size is given, n, the Yule process 
should be conditioned on having n tips: such conditioned branching processes have received sig- 
nificant attention in recent years, due to e.g. Aldous and Popovic [2005], Gernhard [2008], Mooers 
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et al [2012], Stadler [2009, 2011], Stadler and Steel [2012]. This "tree-free" approach for com- 
parative phylogenetics was previously addressed by Sagitov and Bartoszek [2012] and Crawford 
and Suchard [2013], [much earlier Edwards, 1970, used a related branching Brownian process as 
a population genetics model]. 

In our work we show that a properly scaled sample variance is an unbiased and consistent 
estimator of the compound parameter p 2 = o 2 /X which we call the evolutionary rate of the trait 
value in question. Our main mathematical result, Theorem 2.1, gives an asymptotical expression 
for the variance of the phylogenetic sample variance. This result leads to a simple asymptotic 
formula for the estimated standard error of our estimator. Our result is in agreement with the work 
of Crawford and Suchard [2013] whose simulations indicate that their approximate maximum 
likelihood procedure yields an unbiased consistent estimator of o 2 . This is illustrated using the 
example of the Carnivora order studied previously by Crawford and Suchard [2013]. 

The phenotype modelled by a Brownian motion is usually interpreted as the case of neutral 
evolution with random oscillations around the ancestral state. This model was later developed into 
an adaptive evolutionary model based on the Ornstein-Uhlenbeck process by Felsenstein [1988], 
Hansen [1997], Butler and King [2004], Hansen et al. [2008], Bartoszek et al. [2012]. The tree- 
free setting using the Ornstein-Uhlenbeck process was addressed by Bartoszek and Sagitov [2012] 
where for the Yule-Ornstein-Uhlenbeck model, some phylogenetic confidence intervals for the 
optimal trait value were obtained via three limit theorems for the phylogenetic sample mean. Fur- 
thermore, it was shown that the phylogenetic sample variance is an unbiased consistent estimator 
of the stationary variance of the process. 

At the end of their discussion Crawford and Suchard [2013] write that as the the tree of life 
is refined interest in "tree-free" estimation methods may diminish. They however indicate that 
"tree-free" estimates may be useful to calculate starting points for simulation analysis. We cer- 
tainly agree with the second statement but believe that development of "tree-free" methods should 
proceed alongside that of "tree-based" ones. 

One of the most useful features of the tree-free comparative models is that they offer a natu- 
ral method of tree growth allowing for study of theoretical properties of phylogenetic models as 
demonstrated in this work [and also Sagitov and Bartoszek, 2012, Bartoszek and Sagitov, 2012, 
Bartoszek, 2014, Crawford and Suchard, 2013]. Another alternative to studying properties of these 
estimators is the tree growth model proposed by Ane [2008], Ho and Ane [2013], Ane et al. [2014]. 
In this setup the total height of the tree is kept fixed and new tips are added to randomly chosen 
branches. These two approaches seem to be in agreement, at least up to the second moments, since 
e.g. they agree on the lack of consistency of estimating Xq. In Sagitov and Bartoszek [2012] we 
showed that under the Yule Brownian motion model Var [X n ] — > 2a 2 . 

In a practical situation "tree-free" methods can be used for a number of purposes. Firstly as 
pointed out by Crawford and Suchard [2013] they can be useful for calculating starting points for 
further numerical estimation procedures or defining prior distributions in a Bayesian setting. Sec- 
ondly they have to be used in a situation where the tree is actually unknown e.g. when we are 
studying fossil data or trying to make predictive statements about future phenotypes, e.g. develop- 
ment of viruses. Thirdly they can be used for various sanity checks. If they contradict "tree-based" 
results this could indicate that the numerical method fell into a local maximum. 
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The paper has the following structure. Section 2 presents the model, the main results and 
an application. Section 3 states two lemmata and a proposition directly yielding the assertion of 
Theorem 2.1. Proposition 3.1 deals with the covariances between coalescent times for randomly 
chosen pairs of tips from a random Yule n-tree. The properties of the coalescent of a single random 
pair were studied previously by e.g. Steel and McKenzie [2001] and Sagitov and Bartoszek [2012]. 
In Section 4 we state two lemmata needed for the proof of Proposition 3.1. Section 5 contains two 
further lemmata and the proof of Proposition 3.1. In Section 6, 7, and 8 we prove the lemmata 
from Sections 3, 4, and 5. Appendix A contains some useful results concerning harmonic numbers 
of the first and second order. 

2 The main results 

The basic evolutionary model considered in this paper is characterized by four parameters (A , u,Xq, a 2 ) 
and consists of two stochastic components: a random phylogenetic tree defined by parameters 
(A, n) and a trait evolution process along a lineage defined by parameters (Xq, o 2 ) . The first com- 
ponent, species tree connecting n extant species, is modelled by the pure birth Yule process [Yule, 
1924] with the birth (speciation) rate A and conditioned on having n tips [Gernhard, 2008]. For 
the second component we adapt the approach by assuming that for a given i = 1 , . . . , n, the current 

(n) 

trait value X> has evolved from the ancestral state Xo according to the Brownian motion with the 
local variance o 2 . 

(n) (n) 

Treating the collection of the current trait values (X| , . . . ,X„ ) generated by such a process as 
a sample of identically distributed, but dependent, observations, we are interested in the properties 
of the basic summary statistics 

x n = x i H) + - +x!r \ S 2 n = ^-±(X^-X n ) 2 , 

n n i . =1 

the sample mean and sample variance. 

According to [Sagitov and Bartoszek, 2012] we have 

r 9-, ( n + 1 n \ o 2 

see Fig 1, left panel (all simulations are produced using the TreeSim [Stadler, 2009, 2011] and 
mvSLOUCH [Bartoszek et al., 2012] R packages). It follows that the normalized sample variance 

gives an unbiased estimator of the compound parameter p 2 := ^ for the Yule-Brownian-Motion 
model, see Fig 2. In the comparative phylogenetics framework the ratio p 2 can be called the 
evolutionary rate as it measures the speed of change in the trait value when the time scale is such 
that we expect one speciation event per unit of time and per species. The next theorem is the main 
asymptotic result of this paper, illustrated by Fig 1, right panel. 
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Figure 1: Left: True and simulated values of E , right: simulated values of Var \S 2 ] with limit 
equalling n 2 /6 + 1. Each point comes from 10000 simulated Yule trees and Brownian motions on 
top of them. Parameters used in simulations are X = 1, Xq = 0 and o 2 = 1. The grey line on the 
right panel fits a curve based on the convergence rate 0{nT x logn 2 ). 



Theorem 2.1 Consider the sample variance S 2 for the Yule-Brownian-Motion model with param- 
eters (A,h,Xq, (J 2 ). Its variance satisfies the following asymptotic relation 



Tl 2 



Var[S 2 /p 2 ] = 1 + — + 0(n- 1 log 2 n), n ->• ~. 



6 

In terms of our estimator (1), Theorem 2.1 yields 

var ^ /p2 H (log ^;-2 )2 +o( "'''- 

where y = 0.577 is the Euler constant, implying that R 2 is a consistent estimator of the evolutionary 
rate p 2 . It follows that for large n, the standard error (estimated standard deviation) of the unbiased 
estimator R 2 can be approximated by 



bbWJ V + 6 logn+y-2 logn-1.423 ^ {I) 

The estimator of Eq. (1) should be compared to the approximate maximum-likelihood esti- 
mator for the local variance o 2 recently proposed by Crawford and Suchard [2013] in the same 
framework of the Yule-Brownian-Motion model. The main difference between two approaches is 
that in Crawford and Suchard [2013] it is assumed that one knows both the number of tips and the 
total height of the otherwise unknown species tree. The Crawford-Suchard estimator is based on 
a closed form of the distribution of phylogenetic diversity - the sum of branch lengths connecting 
the species in a clade. 
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1L_ 1_ ;L 

n = 1 50 n = 1 50 n = 200 

Figure 2: Histograms of R 2 for left to right top n = 5, 10, 50 and bottom n = 100, 150, 200. Param- 
eters used in simulations are X = 1,X 0 =0 and o 2 = 1. 



As an application of their estimator, Crawford and Suchard [2013] study different families 
of the Carnivora order, estimating o 2 for each of the 12 clades. The data for the log-body- size 
disparities was taken from the PanTHERIA database [Jones et al., 2009]. The data summary and 
the Crawford- Suchard estimates are shown in the left part of Tab. 1 . In the right part of Tab. 1 we 
present our estimates p 2 for the evolutionary rate parameter p 2 = o 2 /X for each of the 12 families 
in the Carnivora order. The standard error is computed using (2). We note that the data does not 
take into account the newly described species Bassaricyon neblina from the Procyonidae family 
[Helgen et al., 2013]. 

In the next-to-last column we list the ratios demonstrating a surprisingly good agreement be- 
tween our and Crawford-Suchard estimates. The ratio is taken between two products: p 2 u n on 
one hand, and & 2 t n on the other. Here u n = E [U n ] is the expected age of the conditioned standard 
Yule process with X = 1, while t n is the clade age assumed to be known in the Crawford-Suchard 
framework. Both p 2 u n and o 2 t n estimate the same quantity - the variance in the trait values for 
the evolution of the corresponding clade. Therefore, one should expect these ratios to be close to 
one. And indeed, the 12 ratios have mean 0.97 and standard deviation 0.20. 

Our estimator and its standard error are computed by simple formulae given above. A major 
weakness of our estimator is relatively big standard error for realistic richness values, see the 7th 
column in Tab. 1 . This can be explained by the fact that we do not use an additional information 
about the species tree, like the height of the tree used in the Crawford-Suchard estimator. 

This close agreement is obtained despite a number of features that complicates the comparison 
between two methods. Our approach in its current form does not allow to take into account the fact 
that some trait values are missing. We calculated p 2 for the trait disparity as if it was computed 
using all n trait values. Moreover, it is not be clear how to take into account the measurement 
variance. As shown by Hansen and Bartoszek [2012] even with a known tree, the measurement 
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Family 


n 


hi 


Disparity 


o yoc.) 




P 


p 2 u„ 
a 2 t„ 


P 2 

& 2 /X 


Felidae 


40 (7) 


33.3 


1.588 


.080 (.009) 


4.279 


.649 (.466) 


1.042 


0.560 


Viverridae 


35 (6) 


37.4 


0.662 


.029 (.004) 


4.147 


.284 (.217) 


1.086 


0.676 


Herpestidae 


33(4) 


25.5 


0.482 


.030 (.003) 


4.089 


.211 (.166) 


1.128 


0.485 


Eupleridae 


8(0) 


25.5 


0.916 


.079 (.010) 


2.718 


.758 (1.72) 


1.023 


0.662 


Hyaenidae 


4(0) 


32.2 


0.805 


.122 (.005) 


2.083 


.999 (19.5) 


0.530 


0.565 


Canidae 


35 (3) 


48.9 


0.678 


.030 (.004) 


4.147 


.290 (.221) 


0.825 


0.667 


Ursidae 


8(0) 


42.6 


0.303 


.024 (.002) 


2.718 


.251 (.569) 


0.667 


0.722 


Otariidae 


16(2) 


24.5 


0.386 


.028 (.003) 


3.381 


.227 (.274) 


1.119 


0.559 


Phocidae 


19(0) 


24.5 


0.751 


.052 (.005) 


3.548 


.410 (.438) 


1.142 


0.544 


Mephitidae 


12(3) 


32.0 


0.570 


.039 (.005) 


3.103 


.384 (.588) 


0.955 


0.679 


Mustelidae 


59 (10) 


27.4 


2.263 


.126 (.014) 


4.663 


.811 (.497) 


1.095 


0.444 


Procyonidae 


14(1) 


27.4 


0.531 


.037 (.004) 


3.252 


.332 (.444) 


1.065 


0.619 



Table 1: Data summary. 2nd column: clade richness (number of missing trait values); 3rd column: 
the clade age in millions of years; 4th column: ■ S 2 trait disparity ; 6th column: the expected 
age u n = E [U n ] of the conditioned standard Yule process with X = 1. 



error can cause very diverse effects. Therefore we would expect the situation to be even more 
interesting when we integrate the phylogeny out. 

In their work Crawford and Suchard [2013] estimated the overall speciation rate to be X = 
0.069 per million years. The last column of Tab. 1 demonstrates that using this common value 
for the speciation rate X produces huge discrepancy between our estimates p 2 for the rates of 
evolution p 2 = o 2 /X and the rates of evolution computed using the Crawford-Suchard estimates 
for o 2 . This observation points out that a fair direct comparison of p 2 and a 2 /X would requires 
specific estimates of the speciation rate X for each of the 12 clades. 

3 Outline of the proof of Theorem 2.1 

We start with a general observation, Lemma 3.1, concerning the sample variance 

I— I 

of n, possibly dependent and not necessarily identically distributed, observations (Y\, . . .,Y n ) with 

n 

sample mean Y = n~ l £ Y,. 

;=1 

Lemma 3.1 If (Wi,W2,W3,W4) is a random sample without replacement from random values 
(Yi,...,Y n ), then 

Var [D 2 ] = Cov [W?, W}] - 2Cov [W 2 , W 2 W 3 ] + Cov [WiW 2 , W 3 W 4 ] + n~ l B n , (3) 
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where 



\B n \ <E [Wf] +4E [W?W 2 ] +E [wfflf] +6E [w}W 2 W 3 ] + 4E[WiW 2 W 3 W 4 ] 
Observe that in terms of the sample variance for the scaled trait values 



Yi := Y, 



(«) 



-X 0 



o/VX 



(4) 



~ 0.2 £j2 

we have S„ = — jp 1 , and to prove Theorem 2.1 we have to verify that 



Varto 2 ] =\ + 7 L + o(n- l \og 2 n). 
6 



(5) 



The Yule n-tree underlying the set of scaled values (4) has unit speciation rate. We call it the stan- 
dard Yule n-tree, and denote by ^ be the a-algebra generated by all the information describing 
this random tree. Under the Brownian motion assumption the trait values (4) are conditionally 
normal with 

E[Y i \& n } = 0, Var\Yi\& n ] = U n , 
where U n is the height of the standard Yule n-tree, see Fig. 3. Moreover, see Section 6, we have 

Cov[^,F;|^] =U n -xjf, (6) 

(n) 

where T^- is the backward time to the most recent common ancestor for a pair of distinct tips (i, j) 
in the standard Yule n-tree, see Fig. 3. For a quadruplet (i, j, k, I) of tips randomly sampled without 
replacement out of n tips in the standard Yule rc-tree, we denote 



(«) 



>) 



Jn) 
L ik ' 



Jn) 



Jn) 



Jn) _ » 



Jn) 



Jn) 
L kl ■ 



(7) 



t ~~ V'* ' 5 ~~ k j7 ' "6 

Lemma 3.2 Le? (Wi , V^2) W3, W4) be a random sample without replacement of four trait values out 
of n random values defined by (4). Then in terms of the coalescent times (7) we have 

Cov [Wf ,W 2 2 ] -2Cov \W},W 2 W 3 ] + Cov[WiW 2l W 3 W 4 ] 



= 2Var 



4 Cov 



Jn) Jn) 



+ 3 Cov 



(»)' 



In view of Lemmata 3.1 and 3.2 which are proven in Section 6, to verify (5) it suffices to show the 
following asymptotic result. 

Proposition 3.1 Consider the coalescent times (7). As n — > oo ; 

P 2 



Var 



'1 



Cov 
Cov 



Jn) Jn)' 
1 ' l 2 

tj , 1.3 



7T 

"6 

2- 
3 



+ 0(n _1 log 2 n) 



7T 
~6 
5^ 2 



18 



+ 0(n- l \og 2 n), 
+ 0(n~ l log 2 n). 
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v. 




0 



2 



3 



4 



5 



Figure 3: An example of a standard Yule n-tree with n = 5. The tree height is U n = T\ + . . . + T n , 
where 7} are the times between the consecutive speciation events. The 10 pairwise coalescent 
times tjy for the tips of the tree are = T5, T$ = T5 + T4, z[ 5 J = %y) = T5 + T4 + T3, and 



Notice that the key Proposition 3.1 concerns only the first component of the evolutionary model 
we study - the standard Yule rc-tree. For the standard Yule n-tree it is well known that the times be- 
tween the consecutive speciation events (7\,. .. , T n ) are independent exponentials with parameters 
(1, . . . ,«) respectively, see Fig. 3. As shown in Gernhard [2008], this property corresponds to the 
unit rate Yule process conditioned on having n tips at the moment of observation, assuming that 
the time to the origin has the improper uniform prior [see also Feller, 1971]. 

4 Coalescent indices of the standard Yule w-tree 

Following the standard Yule n-tree from its root toward the tips we label the consecutive splittings 
by indices 1 , . . . , n — 1 : splitting k is the vertex when k—l branches turn into k branches. We define 
three random splitting indices (as we interested in four randomly chosen tips out of n available): 

• K n is the index of the splitting where two randomly chosen tips coalesce, 

• L n be the index of the splitting where the first coalescent among three randomly chosen tips 
takes place, 

• M n be the index of the splitting where the first coalescent among four randomly chosen tips 
takes place. 

To avoid multilevel indices in the forthcoming formulae, we will often use the following notational 
convention 

KL n : = Kl u , LM n : = Lm„ , KLM n : = Klm„ ■ 

To illustrate these indices, turn to the Fig. 3. If the two randomly chosen tips are (1,2), then 
K n = 4. If the three randomly chosen tips are (2,3,4), then L n = 4, Ki n = 2. If the four randomly 
chosen tips are (2,3,4,5), then M„ = 3, Lm„ = 2, Klm„ = 1- 




35 — 



T 5 + T 4 + T 3 + T 2 . 
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The importance of these random indices comes from the following representations. Denote 

(n) 

U k := Tk + \ + . . . + T n the sum of adjacent times between splittings in the Yule tree. Clearly, 

Jn) d_ Jn) Jn) Jn) d_ (n) Jn) Jn) d_ (n) Jn) Jn) £ Jn) R 
T l — U K„ ' T l AT 2 L„ ' T l v T 2 KL n i % \ v T 3 KLM n ' V°J 

To prove Proposition 3.1 we need to know the distributions of these random splitting indices. The 
next two lemmata giving these distributions are proved in Section 7. 

Lemma 4.1 Then 



P(K n = k) 



n+l 



l,...,n-l. 



,n- 1, 



(k+l)(k + 2)' 
P(L _ h) _ (n+l)(n + 2) 6(k-l) 

(„-l)(„-2) (*+l)(* + 2)(* + 3)' ' 
PfM _ /:) _ (n+l)(n + 2)(n + 3) 12(fc-l)(fc-2) 

/cj (n-l)(n-2)(n-3) (* + 1)(* + 2)(*+ 3)(* + 4) ' 

Lemma 4.2 The random numbers KL n ,LM n ,Kjj^ n have the following distributions 

n+1 12 ( H -±l-l\ k=l n-2 

(n-l)(n-2) (ifc+l)(Jfc + 2)U + 3 /' ' 

(n + l)(n + 2) 72(*-l) /n + 3 





= k) 




= k) 


P ( K LM„ 


= k) 



(n-l)(n-2)(n-3) " (*+l)(Jfc + 2)(* + 3) U + 4 /' * 2 ' 
n + 3 72 /(n + lVw + 2) \ 



,n-2, 



" + 2 144 / n+l \ 

(n-l)(n-2)(n-3)"(*+l)(* + 2)U + 3 /' 



5 Proof of Proposition 3.1 

In view of (8), the harmonic numbers 



(9) 



k=l 



k=\ 



play an important role in our calculations as 



U, 



(n)\2 



H n — H k + (H n — H k ) 2 =H n — 2H n H k + H n + H k —H k 
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Lemma 5.1 We have 



E 
E 



.(«)' 



(n)\2 



>)>)' 



T (»)>)' 



H n -E[H Kn ] , 

H 2 - 2H n E [JfrJ + + E [H 2 Kn - H Kn ] 



H n — H n E 



H n — H n E 
1 



2H Ln +4H KLn 



+ H n + E 



3H m „+5H L m„ + 10H klm „ 



2H Ln H KLn +H K t - 1H Ln - H KLn 



+ H n --E[H Mn ] 



+ g E [Hm„Hlm„ + 3-HM n liKLM n + 4HLM n Hf(LM„ + 2H KLMn — AHlm„ — 2Hklm„] ■ 

In view of Lemma 5.1, proven in Section 7, the asymptotic results stated in Proposition 3.1 are 
computed using the following relations involving the harmonic numbers (9). 



Lemma 5.2 We have as n — > °° 



11 



E[H Kn ]=2 + 0(n- l \ogn) 1 E [H Ln ] = 3 + O^ 1 logn), E [H Mfi ] = — + O^ 1 logn), 
E [H KLn ] = 3 - + 0{n- 1 logn), E [H LM „] = 7 - + O^ 1 logn), E [H KLMn ] = - + 0^ logn). 



Lemma 5.3 Let a n 



a stand for a n = a + 0(n 1 log 2 n) as n — )■ °°. TTzen 

E [Hm„] 



E[^ n ]=^-2, E[#J- 3 



ft 2 , 211 
3 "r 18 ' 

it. _ M 
3 18' 



and 



and 



E[%J: 



► _ 1 

> 2 4' 
jfi _ 15 
2 4 ' 



E [#LM„] : 
39 7t 2 



> 167 n: 2 
'18 3 ' 
85 _ 7^ 
18 3 ' 



E [^KZAfJ 

E [#klm„] : 



, 27T 2 41 
^3 9 ' 

27T 2 49 

3 9 ' 



E [H Ln H KLn ] =4 ^ - E [H M „H LM „} =4^-3-' 

E [Hm„Hklm„] =4 3 — E [Hlm„Hklm„] =4 ^9- — 3-. 



The proofs of the last two lemmata are given in Section 8 using the auxiliary results from Appendix 
A. 

With Lemmata 5.1 - 5.3 at hand, the remaining proof of Proposition 3.1 is straightforward. The 
first statement 



Var 



-IE 



H n + E[H 2 Kn -H Kn ]-(E[H Kn ]) : 



n- 



6 ' 
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is obtained applying the classical relation H n = ^- + 0(n ). Further, Lemma 5.1 yields 



Cov 



Jn) Jn) 



An) Jn) 
1 L 2 



.(«)' 



//„E 



2^„ - 



2H u +4H KLn 



+ H n + -E [2H Ln H KLn + H 2 KLn - 2H Ln - H KLn ] - (E [ZfrJ ) s 



where according to Lemma 5.2 



2H. 



K„ 



2H Ln +4H KLn 



0(n ^ogn). 



Thus, applying Lemma 5.3 we obtain the second statement 



Cov 



>) M 
l ' ^ 



7T 2 1 



6 + 3l + 2 V4 2 



39 7T 2 \ 7T 2 9 3 Tt 1 15\ tt~ 

i i i 



+ 



2 4 



Finally, the third statement follows from 



Cov 



Jn) Jn) 



//„E 



2//, 



3//m„ + 5Hlm„ + 10Hklm„ 



+ H n -(E[H Kn ]) 2 --E[H Mn ] 



+ - E [HM„H LMn + 2H Mn H KLM „ + 4H LM „H KLMn + 2H\ 1Mn — 4-H LMn — 2H KLM „] 



Indeed, according to Lemma 5.2 



2H Kn 



3H m „ + 5H L m„ + 10H klm „ 



0(n ^ogn). 



Moreover, from the following three limits 



1 



H n -(E[H Kn ]) 2 -^E[H Mn ] 
E [Hm„Hlm„ + 2Hm„Hklm„ + 4-Hlm„Hklm„] 
E [ h klm„ ~ 2H LMn - H K LM„] 
we get the stated overall limit 

% 2 185 1 /1349 137T 2 \ 2/27T 2 77 



K 2 185 
18 ~^4"' 
1349 137T 2 



18 54 + 9V 18 



+ 



18 3 
27T 2 77 



5n z 



9\ 3 



18 
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6 Proofs of Lemmata 3.1 - 3.2 

Proof of Lemma 3.1. Using the representation 
D 



we find that 



EM=^(lE[^] + ££E[F/y,?]) 



- ^huiZZvm +I£e [Y iY f] +EE E e ) 
+ ;? t^(EEe +EE E E [ y iVt] +EE E 4^1 

If (V^i , W2, W3, W4) is a random sample without replacement of four out of n trait values, then 

EK]=„-'£E[tf], 



e [^] = ^T)EEe[^.] 



e [^] = ^ttEEe[^/] 



and 



E I W > W l = ,(,-!)(. -2) EE E e [n V*] • 

E[ Wl W 2 W 3 W 4 ) = n{n _' 2)(n _ 3) i:i: E I E [F ;W/ ] 



Therefore, we have 



E M = n^E -^E [Wf W 2 ] + "~T~~TV~E [W?W 2 2 ] 

— 1 J 

- E [Wf^] + ("-2)(n-3) E [WlW3W4] . (10) 

n(n— \) n[n—l) 
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Since 



E[D 2 n ] = E[W?]-E[W l W 2 ] 



we conclude 



Var[D 2 ] = n _1 E[Wf| -4n _1 E [W?W 2 ] +Cov [wf,W 2 2 ] - -jj- — jjEfWfW, 2 ] 



- 2Cov \W?, W 2 W 3 ] + 2 ^ 3n 4 ) E \W?W 2 W3] 
L J n[n— 1) 

+ Cov [WiW 2 , W 3 W 4 } - 2 ^ ~^ E [Wi W 2 W 3 W 4 } 



The stated relations follow with 

B n =B[w}] -4E [W?W 2 ] 
2w2 



^— ^ E [wfwfl + 2< " 3n E [Wf W2W3I - 2< ^ 2W 7^ E [Wi W2W3W4] 
« — 1 « — 1 .,i 



□ 



(n 

Proof of Lemma 3.2. Denote by Y-< ' the normalized trait value of the most recent common 



ij 



ancestor of the tips (i, j). Let stand for the a-algebra generated by the pair Yfj 1 '), then 



E 
Var 
Cov 



(«) 



(«)■ 



y(«) 

0' ' 



(«) 



Var 



(«)' 



3>W 



(„) 



0, 



implying (6) 



'./ 



Cov [F;,F^] = Var 
By Eq. (13) of Bohrnstedt and Goldberger [1969], we have 

Cov [Z;Zy, Z k Zi] = mim k c jt + mmc^ + m jm k cu + m jmic ik + c ik c jt + c a c jk 

for any sequence of normally distributed random values Zi,Z 2 ,... with means E[Z;] = m ; - and 
covariances Cov [Z,-,Z 7 ] = c,j. In the special case with m ; - = 0 it follows 

Cov [ZiZj,Z k Zi] = c ik Cji + cuc jk , 
Cov [zf,ZjZ k ] = 2cijCi k , 
Cov [Zf,zf\ =24. 
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(n) 

Using conditional normality of Y, and putting c, 7 = U„ — T>- ; , we derive from these relations that 



(»)\2 



Cov[Y^Yf\%]=2(U n -x};>) 



Cov [Y?,YjY k \9r n ] = 2(U n - xf){U n - zf ), 



Cov [7^,7,7;!^] = ([/„ - Tf)([/ n - T W) + (f/ n - T W)([/ n _ 
yielding in terms of (7), 

Cov [W?,Wf|^] = 2(t/ n -T 1 (n) ) 2 , 



COV [Wf , W 2 W 3 |%] = 2(U n - Tr^f/n - X 



2 /' 



Cov [WiW 2 ,W 3 W 4 \y n } = (U„ - 4 n, )(u n - xf) + (U„ - x^){U n - T 4 
By the total covariance formula, we derive 



Cov [wf , W 2 2 ] = 2E \(U n - x 
Cov [W?, W 2 W 3 ] = 2E [(t/ n - x 
Cov [Wi W 2 , = 2E [(£/„ - T 



3E 



(U n -x 



n )\2 



+ Var[t/„] 



] )(U n -x 



+ Cov 
+ Cov 

E 



U TJ — x y 



u„ 



x {n) U 



Combining these relations we get 

Cov [Wf , ] - 2Cov [wf , W2W3] + Cov [Wi W 2 , W 3 W 4 ] 

'(U n -x[ n) )(U n -x^ 



2E 



+ 3E 
+ Var[t/„]-2Cov 



4E 



(C/n-rf n) )(C/n-^ n) ; 



r/ u — t v 



u„ 



This together with 

Var[t/„]-2Cov \u n ,U n -x[ n) 



- E 



(«) 



= E[t/, 2 ]-E[C/„] 2 -2E U n (U n -x\ n) ) +2E[f/„]E 



E[t/ 2 ] -2E U n (y n -x\ 



(«) 



(„) 



= E 



(»)\2" 



«J\2 
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implies the assertion of the Lemma 3.2 

Cov \W?,W%] -2Cov [Wf ,W 2 W 3 ] +Cov[WiW 2 ,W 3 W4] 



= E 

-4E 
= E 

-4E 
= 2Var 



{n)\2 



+ E 



n )\2 



K 1 ) 



(U n -T[ n) )(U n -T { 2 



(»)> 



-E 
+ 3E 



(«)\2 



(U„-T 1 , 

(«)» 



+ Var 
+ 3E 
4 Cov 



,(«) 



+ E 



U n x\ 



(n)l 



+ 3 Cov 



□ 



7 Proofs of Lemmata 4.1, 4.2, and 5.1 

Proof of Lemma 4. 1 From the definition of K n it is easy to see that, for k = 2, . . . , n, 

1 (k+l)(k-2) 



P(K n <k-\\K n <k) = \- 

r 

Therefore, 



P{K n <k-\) 



P{K n = k-\) 



g) *(*-!) 

(n + l)(n-2) n(n-3) (n-l)(n-4) (*+l)(fc-2) 
»(n-l) (n-\)(n-2) (n — 2)(/i — 3) jfc(jfc-l) 

(n + lKfc-2) 
(n-l)fc ' 

(n+l)(Jfc-l) (n+l)(*-2) « + l 2 



(n-l)(fc+l) (n-l)A: /i-l(*+l)fc 
Similarly, for k = 3, . . . , n, 

3 (*+2)(*-3) 



P(L n <k-\\L n <k) = 1 
P(L n <£-l) 



© *(*-!) ' 
» + 2)(n-3) (n + l)(n-4) n(n-5) (* + 2)(*-3) 



P(L„ = £-1) 



n(n-l) (n-l)(n-2) (w — 2)(w — 3) fc(Jfc-l) 
(K + 2)(K + l)(fc-2)(fc-3) 
(tj — 1) (/i — 2) l)ik ' 
(7» + 2)(n + l)(Jfe— l)(ik — 2) (n + 2)(n+l)(Jk-2)(Jfc-3) 



(n-l)(n-2)(fc + 2)(fc+l) (n— — 2)(Jfc+ l)fc 
6(n + 2)(n+l)(*-2) 



(n-l)(n-2)(Jfc + 2)(*+l)*' 
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and, for k = 4, . . . , n, 

6 (* + 3)(Jfc-4) 



P(M n <£-l|M n <£) = l 

P(M„<*-1) 
P(M n = k-l) 



© *(*-!) ' 
(n + 3)(n + 2)(n + l)(fc-2)(fc-3)(jfc-4) 
(/i-l)(w-2)(/z-3)(]fc + 2)(]fc+l)]fc 
12(/i+3)(n + 2)(n+l)(Jfc-2)(*-3) 



(n - 1 ) (n - 2) (n - 3) (Jfc + 3) (k + 2) (Jfc + l)k 



□ 



Proof of Lemma 4.2 Clearly, 



P(K Ln =k)= £ P(K, = k)P(L n = l) 

l=k+\ 

(n+l)(n + 2) 12 "f, 1 1 



(n - 1) (n - 2) (* + 1) (* + 2) ^ (/ + 2) (/ + 3) 
leads to the first assertion. Further, 

(n + l)(7i + 2)(n + 3) 



P(Lm„ = k) 



(n-l)(n-2)(n-3) 

12(m-l)(ra-2)(ra+l)(ra + 2) 6(1-1) 
X m ^ +1 (m+l)(m + 2)(ro + 3)(m + 4)(ro-l)(ro-2) ' l)(/c + 2)(/c + 3) ; 



and 



p / 7 .x (n + l)(n + 2)(n + 3) 72(*-l) "f, 1 

P(L M „ = k) 



(„-l)(n-2)(n-3) (* + l)(* + 2)(* + 3) ^ 1 (m + 3)(m + 4) 
(n + l)(n + 2)(n + 3) 72(*-l) , 1 1 , 



(n-l)(n-2)(n-3) (k+ l)(/c + 2)(/c + 3) K k + A n + 3' 



Finally, 



P(K =k)= 'y (n + l)(w + 2)(n + 3) 72(m-l) m+1 2 

1 LM " 1 M ^ +1 (n-l)(n-2)(n-3) (m+ l)(m + 2)(m + 3)(m + 4) m- 1 (*+l)(* + 2) 

^ 2 (n+ l)(n + 2) 72(m-l) m+1 2 



m i^. 1 (n-l)(n-2)(n-3) (m+ l)(m + 2)(m + 3) m- 1 (*+l)(* + 2) 

(n+l)(n + 2)(n + 3) 144 1 

(n-l)(n-2)(n-3) ' (* + l)(* + 2) m ^ +1 (m + 2)(m + 3)(m + 4) 

(n+l)(n + 2) 144 1 

" („- 1)(„ - 2)(n- 3) ' (fc+l)(fc + 2) m L +l ( m + 2)(m+3)' 
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and therefore, it remains to use the equalities 



n-2 



1 



1 



1 



n-2 



m i r+1 (m + 2)(m + 3) 
2 



k+3 n+l 
1 



w ^ +1 (m + 2)(m + 3)(m + 4) (k + 3)(k + A) (n + l)(n + 2) 



Proof of Lemma 5 . 1 . We have 



□ 



(n) 



5. 



and 



(tD 2 =E (t/W) 2 = J f7„ + // 2 -E[^+2//^-// 2 J 



Further, using (8) and t/|? = t/| n) + U&>, we get 



Thus 



Jn)Jn) 



3 
E 



(^; } ) 2 



4 

+ 3 E 



+ 3 E 



(4t } ) 2 



>)>)' 

1 l 2 



: ff 2 - 2# n E [H Ln ] + H n + E [H 2 Ln -H Ln ]+-E [(H n - H Ln ) (H Ln - H KLn )\ 
+ \ E K - 2 ^^„ + + H 2 KLn - H KLn ] 



H n — H n E 



2H Ln +4H KLn 



+ H n + E 



2H Ln H KLn + Hh - 2H Ln - H KLn 



Finally, for two pairs of sampled tips, we have three coalescent events to consider: going from four 
to three selected nodes, 4 — > 3, going from three to two selected nodes, 3—7-2, and going from two 
to one selected nodes, 2 — > 1. The coalescent 4 — > 3 holds across the two pairs with probability 
7^r = I and within a pair with probability i. Given the former outcome, the coalescent 3—7-2 

holds again across the pairs with probability | and within a pair with probability |. Otherwise, 
the coalescent 3—7-2 holds across the pairs with probability | and within the second pair with 
probability ^. The four possibilities (| x |, | x | , | x |, | x |) produce the following four terms 
in 



>)>) 
n l 3 



/rr(«) \2 
\ KLM„' 



4 

+ 9 E 



LM„ KLM„ 



2 

+ 9 E 



J A*) J A") 

U M„ KLM„ 



1 

+ 5 E 



' T M) TT {n) 
u M n U LM„ 
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It follows, 



>)>) 



2 ( M „) \2±—TT( n hr( M ") I ^ , ^ („) (M„) 

\ U M„ ) g \ KLM n ) q U M„ U ^LM„ q, U LM„ cy /fLM„ q, Ly M„ U LM„ 



Using the representation for E 
1 



r {M n ) jM n ) 
'1 L 2 



= E 



(M„) x 2 
KLM„> 



TJ {M n ) jj(M„) 
KLM„ LM„ 



H M —H, 



2H LMn +AH 1 



KLM„ 



+ H Mn + 



2H LM „H KLMn +H% LM)i - 2H LM „ -H KLM „ 



we can write 



Hn - E [2H n H M „ + H M „ - Hm h ] +H n + -E 



Ha* —H, 



2H LMn +AH l 



M„ - n M» " 



KLM„ 



2 

+ 3 E 



2H LM „H KLMn + H\ LMn - 2H LM „ - HfCLMn 



+ -E[(H n - H M „ ) (H M „ - H LM „ )] + — E [ (H n - H M „ ) (H M „ - H KLM „ ) ] 
which after a rearrangement gives the last statement. 



8 Proof of Lemmata 5.2 - 5.3 



□ 



In this section we will often use the elementary relations of the following type 

6 12 1 

+ ■ 



(k+l)(k + 2)(k + 3)(k + 4) (k+l)(k + 2) (k + 2)(k + 3) (k + 3)(k + 4)' 

(k-\)(k-2) 1 5 5 

(k+l)(k + 2)(k + 3)(k + 4) ~ (k+l)(k + 2) ~ (k + 2)(k + 3) + (k + 3)(k + 4)' 

6k 15 4 

+ 



(k+\)(k + 2)(k + 3)(k + 4) (k+\)(k + 2) (k + 2)(k + 3) (jfc+3)(*+4)' 

k(k-5) 1 6 _ 6 

(k+l)(k + 2)(k + 3)(k + 4) ~ (k+l)(k + 2) + (k + 3)(k + 4) ~ (k + 2)(k + 3Y 

valid for all A: > 1 . 



(ID 
(12) 
(13) 
(14) 



PROOF of Lemma 5.2. The first three stated relations are obtained using Lemmata 4.1 and A.l. 
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Equalities 



r rc+1'^ 1 2H k 2(n-H n ) 

[ " J n- 1^(^+1)^ + 2) n-l 

(n+l)(w + 2) "f. 1 6(*-l)H» 



(n - 1) (n - 2) ^ + !) + 2 ) + 3) 

6(n+l)(n + 2) ^ 2//, fl, ^ _ 3n(n + 1) - 6nfl„ 



(n-l)(n-2) V(fc + 2)(fc + 3) (fc+l)(fc + 2)V (n-l)(n-2) ' 

give the first and the second stated relations, and the third one follows from 

(n+ l)(n + 2)(n + 3) ' y l \2{k-l)(k-2)H k 
L M «J ( n -l)(n-2)(n-3)^ (*+l)(Jfc+2)(jfe + 3)(Jk + 4) 

(12) 12(k + 1)(/z + 2)(k + 3) n y l / fl, 5fl, ggfc 

(n-l)(n-2)(n-3) ||l(Hl)(t+2) (fc + 2)(fc + 3) (fc + 3)(fc + 4) 

__ ^-n 3 + 30n 2 + 5|i-n-20-12(n 2 + l)//„ 
= (n-l)(n-2)(n-3) ' 

The second three stated relations are obtained similarly using Lemmata 4.2 and A.l. Indeed, 

(n+lKn + 2) ^ 1 12//, n+1 "p, 1 12//, 

E[fl*zJ- t- L 



(/»-l)(n-2) J ^ 1 (*+l)(*+2)(* + 3) (n-l)(n-2) fe ti (*+l)(* + 2) 
6(n + !)(« + 2)^ fl]t "f, 1 ^ > 12(n-fl„) 



(„-l)(„_2) V^(*+l)(* + 2) k ^(k + 2)(k + 3)J (n-\)(n-2) 
6n(n-H n ) 6(3n 2 + 5n - 4(n + l)fl„) _ 6fl n 3n(n-5) 



(n-l)(n-2) 4(n-l)(n-2) (n-l)(n-2) 2(n-l)(n-2) 

implying E [Z/axJ = § + 0(n _1 logn). Furthermore, 

err/ l V 1 12{k-\)H k . (11X13) 7 , _ x . 

E fl/Af = / t n — w — r + 0(n logn) = - + 0(« logn ), 

££(fc + l)(* + 2)(Jfc + 3)(* + 4) V S ; 3 V 5 ; ' 

oo V2^^ 4 



□ 



Proof of Lemma 5.3. The stated relations are obtained using Lemmas 4.1, 4.2, A.l, A.2. Firstly, 

E [H y = 5±i y 2 "- 2 = ?&±1> ( a „ + "- tf -- 2 "" 

L " J «-l,=i (*+l)(* + 2) B-l I" «+l 
2(fl„(n + l)+n-// 2 -2fl„) 7T 2 



n-l 3 



+ 2. 
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Similarly, we have 

_ n + 1 1, 1 2// fe _ 2(n+ 1) n//„ - ft _ 2ft// n - 2n 7T 2 
E ™ = L (jt+l)(it + 2) = Ti-l n+l ~ n-l ^ T 

oo 2/"/ 

Observe that the limit is Lr=i (I+T)(i+2) ■ m * ne same manner we obtain 



^(A:+l)(fc + 2)(fc + 3) £i(*+l)(* + 2) £i(* + 2)(* + 3) 
p rr> 1 . r 12/4 _ y 6/4 y 3 , 6Hk 

[ ^ J £l(*+l)(* + 2)(* + 3) l)(* + 2) £i(*+2)(*+3) 

Using the decomposition (12) we find 

F r„ 2 1 _ f. _m-m-m^_ ^ 211 

L Mj £ 3 (*+l)(* + 2)(* + 3)(* + 4) 3 ^ 18' 

" 12(fc-l)(fc-2)ff fc 31 
1 M " J £ 3 (*+l)(* + 2)(* + 3)(* + 4) 3 18' 

Using the difference between (13) and (1 1) we find 

2 " 12{k-\)Hl _167 ^ 2 

L LM " J ^ (fc+l)(£ + 2)(fc + 3)(& + 4) 18 3' 



E [H LMn ] 

Using (11) we find 



12(k-\)H k 85 7T 2 



f 1 (*+l)(* + 2)(* + 3)(* + 4) 18 3 



°° 72/7 2 ?77- 2 41 

p \lt2 ] . y /z -"fc _ 

L ™" J ^fc+l)(£ + 2)(£ + 3)(£+4) 3 9 ; 

Ff/7 I^V ™ fe 49 

L KiM n \^ (Jfe+l)(Jk + 2)(Jk + 3)(ik + 4) 3 9" 

Since E [Hg n ] = 2( „jf"' ) , we have 

F[w „ i-tV 6(m-\)H m 2(m-H m ) 
1 Ln KLnl ^ ^ 2 (m + 1) (m + 2) (m + 3) m - 1 

f, 12m// m ~ 12// 2 



^ (m+l)(m + 2)(m + 3) ^ (m+ l)(m + 2)(m + 3) 
15 tt 2 9 39 Tt 1 
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where we use the following corollary of Lemma A. 1 

2mH m ^ 3H m ^ H„ 



I ^ 1 (m+l)(m + 2)(m + 3) ^(m + 2)(m + 3) ^(m+l)(m + 2) 4 

Similarly, 

P|H „ i_>v 12(m-\)H m 2(m-H m ) 

^ 2 (m+ l)(m + 2)(m + 3)(m + 4) m— 1 

_ y UAmH m y 144fl% 

~ ^ (m+l)(m + 2)(m + 3)(m + 4) ~ ^ (m+ l)(m + 2)(m + 3)(m + 4) ' 



where 



f W 24(-l + 15/4 - 22/9) = 

~ 144//2 (ii) 4tt 2 82 



-(fc+l)(fc + 2)(fc + 3)(fc + 4) 3 9 



so that E [H LMn H KLMn ] -» 9 , . 

Further, in view of E [ff/J = , the limit for E [//m„#lm„] can be computed as 



E[H Mn H LM „] 



(n-l)(n-2) 

^ 12(m— l)(m — 2)// m 3m(m+l)—6mH n 
^ (m+l)(m + 2)(m + 3)(m + 4) (m- l)(m-2) 

P 36mH m ^ 72mH^ 



^ (m + 2)(m + 3)(m + 4) ^ (m+ l)(m + 2)(m + 3)(m + 4) 

where 

^ (m+ l)(m + 2)(m + 3)(m + 4) ~ 36 + 216' 

V~< WlH m r-i 2H m r-i 17 

^ (m + 2)(m + 3)(m + 4) ~ ^ (m + 3)(m + 4) (m + 2)(m + 3) ~ 36' 

221 ^ 2 



yielding E [H Mn H LMn ] =4 ^jf - \. Finally, from 



6//„, 3m(m-5) 
E L-n*zJ = 7- — ^ + ; 



(m— l)(m — 2) 2(m— l)(m — 2) 
we get 

p r„ „ 1 y \2{m-\){m-2)H m / 6H m 3m(m-5) 

L m„ klm„J =4^ (m+l)(m + 2)(m + 3)(m + 4) V(m-l)(m-2) 2(ro- l)(ro-2) 

_ y 72i/£ y 18m(m-5)// m 

~ ~! (m+l)(m + 2)(ra + 3)(m + 4) ^ (m+ l)(m + 2)(m + 3)(m + 4) 
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where 



E 



12H m 



(ii) 2n 2 41 



m=l 



(m+ l)(ra + 2)(ra+3)(ra + 4) 
m(m - 5)i/ m 



3 9 ' 



(14) 1 



I 



(m+ l)(m + 2)(m + 3)(m + 4) 



6' 



m=l 



so thatE[// M „^LM„] =4 



14 
9 • 



□ 
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A Auxiliary results involving harmonic numbers 

Some of the following results can be found in Adamchik [1997] and Sofo [201 1, 2012, 2013]. 
Lemma A.l We have 
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In particular, 
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and generally for m > 1, 
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Using the classical relation 2£™ =1 -g = H^+H m which follows from 



m tj m i k i m i m i m tj tj m tj 

\^ n k V V V V V n m~ n i -1 tj2 , fj \-< n k 

Lt= LrLy-LyLi-l- ; ^+«™ -It- 

fe=l " Jt=l * (=1 4 i=l l k=i K i=l 1 k=l * 



we get 



w v z i + m / n m 



(=1 



V / H n +k-l — H n -\ H n+ k-\—H n i 

2j I 7. 1 j ~ Hnytin+m-l ~ tl n -\) 



k=2 

TJ in TJ tj m l rr rr 
t> i r/2 , r> V- "n+k-] ~ ™n-\ tl n+k -tl n 
— H n +H m -\-H m ^ 7 2^ I- tinytin+m-l — tin-\ , 

m fe=1 /c fe=1 /c 

Thus 

V 1 ^fe ^ / l> i l/2 , f> H m \ H n -i(H n+m — H n ) H n -\H n 

2^ 77~~ \77~~ r~T\ — ~~ I "n + "m + i 

rZi (/c + mj(/c + m+ 1) m \ m/ m n + m 



k=l 

\ m TJ TJ T m ~ 1 ZJ ZJ TJ I TJ TJ 



m k m k m 

To finish the proof it remains to observe that - B± j — - — > 0 as n — > °° for any fixed k. 
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